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Abstract. We introduce a new class of finite differences schemes to 
approximate one dimensional dissipative semilinear hyperbolic systems with 
a BGK structure. Using precise analytical time-decay estimates of the local 
truncation error, it is possible to design schemes, based on the standard 
upwind approximation, which are increasingly accurate for large times when 
approximating small perturbations of constant asymptotic states. Numerical 
tests show their better performances with respect to those of other schemes. 



1. Introduction 

Consider a BGK system in one space dimension, for the unknowns /' G 
k > 1 and i = 1, m: 

f 'A/' + Xid x f = Mi(u) - f\ 

(1) m - 

[ where u := YhLiF- 

Here x € R and t > 0, the A;, for i = 1, ...,m, are given distinct real values, and 
the functions Mi = Mi{u) G M. k are smooth functions of u such that: 

m 

y^ y Mj(u) = u. 

i=l 

Following [3J, we rewrite system (JTJ) in its conservative- dissipative form. This means 
that we assume that there exists an invertible matrix 

Dn D 12 
D 21 D 22 

such that, setting mi = k, m 2 = k(m — 1), the new unknown 
(3) Z = Df = (u, Z) T G R mi x M™ 2 , 

solves the system 

(4) 



(2) D = 



d t u + A u d x u + A 12 d x Z = 0, 
d t Z + A 21 d x u + A 22 d x Z = Q(u) - Z, 



where A is symmetric and Q(u) is quadratic in u, i.e.: Q(0) = and Q'(0) = 
0. Observe that, after the transformation, the source term is zero in the first 
component and the second one is the sum of a quadratic term and of the dissipative 
term — Z. 

It is proved in [13] and [3] that, under some additional conditions, usually induced 
by suitable entropy functions, and for initial data which arc small perturbations of 
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constant equilibrium states and smooth in some suitable norms, the corresponding 
smooth solutions exist globally and their i°°-norm decay, for large times, as 

(5) u = o(t- 1 / 2 ), z^oir 1 ), 

and similar estimates are available for their space and time derivatives. Notice 
that the improved estimate for the unknown Z can only be obtained in these new 
coordinates and does not hold for other combinations of the unknowns. 

The aim of this paper it to take advantage by these time decay estimates to 
build up more accurate numerical schemes. To be more specific, we show in the 
following that for standard numerical schemes, for instance upwind schemes with 
the source term approximated pointwise by the standard Euler scheme, the local 
truncation error for the conservative-dissipative unknowns (u, Z) has the following 
decay as t — > +00, for a fixed CFL ratio: 

7 u {x,t) = 0{Ax t~ 3 / 2 ), 7 2 (x,t) = 0{Ax t™ 3 / 2 ). 

It can be seen numerically that the corresponding absolute errors, for a fixed space 
step, decays as 

e u {t)=0(t- 1 ' 2 ) 1 ez{t) = 0{r 1 ), 

which implies, taking into account ([5J), that the relative error is essentially constant 
in time. 

Here, our main goal is to improve the decay rate of the truncation error to achieve 
an effective decay in time of the relative error, both in u and Z . Before presenting 
our strategy and our main results, let us review some different attempts to design 
effective numerical approximations for hyperbolic equations with a source term. 
Let us mention some families of schemes, sometimes overlapping: Well Balanced 
[H EE [aHUUll [5], Runge-Kutta IMEX HQ] , upwinding source [HJ [2j 1 EE] , and 
asymptotic preserving |14[ 116] . The main idea in all these schemes is to use some 
knowledge of the actual time behavior of the solutions to improve their numerical 
approximation, at least in some specific regimes. 

In particular, in [T], the linear version of the present problem was considered and 
therefore, to approximate the solutions around non constant asymptotic states, 
some schemes were proposed, which had the property to become higher order 
(in space) for large times, thanks to the careful consideration of the analytical 
decay rates of the solutions. A different attempt was given by the well balanced 
schemes, see for instance [121 El E], namely schemes which are exact when computed 
on stationary solutions of the problem, even if up to now, the time decay rate 
of the unsteady solutions has not yet been explicitely considered. However, the 
Asymptotic Preserving properties of some Well Balanced schemes, can yield nice 
results for large times, as in [TT], see Section [7] below. 

In the present striking difference with these previous works lies in the fact 

that the asymptotic equilibrium states are constant and therefore all the consistent 
schemes are exact on them. So, the goal of our work is a bit different. In this paper, 
we design schemes which are able to improve their performance for large times, when 
the initial data are small perturbations of a given constant equilibrium state. To 
obtain these results, we use the estimates in [3] to perform a detailed analysis of 
the behavior of the truncation error for a general class of schemes, which generalize 
and improve those introduced in pQ. Thanks to this analysis, we are able construct 
schemes such that the truncation order behaves as 

7 u {x, t) = 0{Ax t~ 2 ), 7z(x, t)) = 0(Ax t~ 2 ), 

for a fixed CFL ratio and such that their asymptotic numerical error, observed in 
the practical tests, improves of t^ 1 ! 2 on the previous schemes. 
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The plan of the paper is the following. In Section [2J we introduce our analytical 
framework. The main schemes are derived in Section [3J where we show how to 
improve the time decay of their local truncation error. Sections [4] and [5] are devoted 
to the monotonicity conditions for the new scheme in the 2x2 and 3x3 cases 
respectively. Then we present some remarks in the linear 2x2 case, to allow a direct 
comparison with other schemes. Scction[7]prcscnts some numerical tests which show 
the nice behavior of our new schemes both in the linear and the nonlinear cases. 

2. The analytical framework 

Let us observe that when transforming system ([T]) in system we can always 
assume that the block Dn and D12 have the special form 

D 1X = I k , D u = (J fc J fc • • • I k ) e R fcx ™ 2 , 

and setting A = diag{\\Ik, • A m /fc), we have that 

Therefore, we can rewrite our system in a more compact form: 

(6) d t Z + Ad x Z = -Z + DM(u). 

To guarantee the existence of the matrix D in $2$, we can assume that our 
system is strictly entropy dissipative in the sense of [13] and verifies the Shizuta- 
Kawashima condition [HI [T3J |3J. For instance, following Bouchut [J], we may 
assume the following sufficient condition. 

Condition ED (Entropy Dissipation condition) There exists an open set 
il C R fe and a strictly convex function r\ = 77(14) : i7 — > R ; such that the 
matrix Ml(u) T 77" (u) is symmetric and strictly positive defined for all u G ft and 

i = 1, ... ,771. 

Under this condition, and using the results in [H [T31 [3] , it is possible to prove 
the existence of a matrix D in (j2Jl, with all the properties mentioned above. The 
details of the derivation of the actual coefficients of this matrix, which in general is 
not unique, but depends on the entropy dissipative function, can be found in the 
general case in [3] , and are not relevant for the following discussion, even if of course 
they are relevant to derive the final schemes. However, in the examples presented 
below, the matrix D is always explicitly given. 

2.1. Algebraic conditions. Let us state now some general properties for the 
matrix D that we shall use in the following, and which are easily derived under 
condition ED. 

First of all, since A = DAD^ 1 in ^ is symmetric, we have that the matrix 
H := D T D and the diagonal matrix A commute, i.e. 

(7) D T DA = AD T D, 

So, for the generic ij block element of the matrix H, we have that (A-s — Xj)Hij = 0, 
hence 

Hij = for i ^ j once A; 7^ Xj . 

Since A is a block diagonal matrix with m distinct eigenvalues with multiplicity 
egual to k, we get that H is a block diagonal matrix of the following form 

(8) H = D T D = diag(hi, h m ), h l eR kxk , i = l,..,m. 



4 



D. AREGBA-DRIOLLET, M. BRIANL AND R. NATALINI 



Moreover, the matrix H = D T D is symmetric and invertible and we can find the 
following expression for D^ 1 , 

(9) D- 1 = H~ l D T = ( H ^ H i~ ljD 2i 



H 2 D 12 H 2 D 



22 



where H\ = hi and H 2 = diag(Ii2, .., h m ). Finally, using (|8|), the following relations 
hold: 



(10) D 12 = -D^D 22 , D^D 21 = J?i - I k , D{ 2 D 12 + D L 22 D 22 = H 2 . 

There is another important relation involving the matrix H, which will be of use in 
the following. Let An be the top-left block of the symmetric matrix A in system 
@. Therefore: A n = AxH^ 1 + D 12 K 2 H 2 X Dl 2 , i.e. 



(11) A 11 ^X 1 h 1 1 + (lkh---Ik) diag{\ 2 h 2 1 ,...,\ rri hn) 



( Ik \ 

h 



V h J 



i=l 



From now on, we shall consider each matrix N £ ^kmxkm ^ Q ^ e decomposed 
into blocks as follows: 



(12) N 



N 2 i N22 



with Nn 6 M fexfe , N 12 € R fexm2 , N 2 i G R m2Xfc ,7V 22 € M m2Xm2 . For the particular 
case of a block diagonal matrix N = diag(ni, ...,n m ), m £ ~R kxk for i = 1, ...,m, 
we shall use the following notation 

(13) -/V = diag(N±, N 2 ), N\ = n\ and A^2 = diag(n 2 , ..,n m ). 
Moreover, for a generic vector V € R fcm we shall write 

(14) V=(Vi,V), with Vl G V G M m2 . 

Example 2.1. Consider the special case of the 2x2 hyperbolic Jin-Xin relaxation 
system (171115]: 



(15) 



d t u + d x v = 0, 

d t v + X 2 d x u = F(u) 



for A > 0, where the unknowns u and v are scalar and the function F = F(u) is 
smooth, with F(0) = 0. This case is obtained from ([T]) for k — 1, m = 2, and 
A2 = — Ai = A, by setting 

u = f 1 + f, v = X{f - J 1 ), F{u) = A(M 2 - Mr). 

Recall that for bounded initial data, and for A > M, where M is a positive 
constant which depends on F and on the initial data, there exists a global bounded 
solution to the Cauchy problem ifTB]) . see [H]. Under the weaker condition 

(16) A>|F'(0)|, 

the problem is dissipative in the sense of [13| and the Shizuta-Kawashima condition 
is verified, at least for small values of u, and so, at least for smooth and small initial 
data, there exists again a global smooth solution to the same problem. To obtain 
the time decay rates of these solutions, we need to rewrite the problem in the 
conservative-dissipative coordinates. Assume Q16p and set a = F'(0) so that the 
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quantity /i = (A 2 — a 2 ) 1 I 2 is real and positive. Setting Z = fj,(v — au), the new 
unknown (u, Z) solves the problem in form Q: 

( d t u + d x (au + ±Z) =0, 

(17) { 

[ 8 t Z + d x {% - aZ) = n{F{u) - au) - Z. 

In this case the matrix D is given by 

1 1 

where a± = A ± a > 0, from assumption (|16[) . 



D = 



Example 2.2. Let us now compute the conservative-dissipativc form for the 
following 3x3 BGK model, 

f d t fi — Xd x fi = Mi{u) - fi, 
(18) ^ d th=M 2 {u)-f 2 , 

{ d t f 3 + Xd x f 3 = M 3 (u)-f 3 . 

Let F = F(u) be a smooth scalar function such that F(Q) = and let 7 be such 
that j'{u) = \F'(u)\, with 7(0) = 0. We choose our three maxwellian functions as 
follows, for (3 e]0, 1[ and A > 



M 



(19) 



(u) = I (iM^M + ^ , m 3(u) = 1 (tM+£M + ^ , 



M 2 (u) = u - Mi(u) - M 3 (u) = (1 - (3)u - 2£2. 



The functions M^, i = 1,2, 3, are strictly increasing if for any u under consideration 

\F'{u)\ 



(20) 



A > 



1-/3 ' 

and so condition ED is verified. Let a = F'(0) and a = \a\ + /3A, following the 
results in [H [TjJl [3] it is possible to compute the matrix D for the change of variables 
© as 

/ 1 1 1 \ 



q+q / A(q — a) 







A(a — a) 
a(a+a) 



A-a 



A-a 



A-a 



2.2. Time Decay Properties. Here we report some time decay results which 
have been mainly proved in [3] , and which will be useful in the following. 

Theorem 2.3. Let Z = (u, Z) be the local smooth solution to the Cauchy problem 
for system @. Let E s = max{||Z(0)||^i, ||Z(0)||ffs} a norm for the initial data, 
which are taken in H k for k large enough, and assume E 2 small enough. Therefore, 
under the Condition ED, the solution is global in time and the following decay 
estimate holds, for all j3 < k: 

(21) \\d^Z{t)\\ L ^ <Cwm{l,t-*-Z}Ep +l , 

with C = C(Ep +a ) for a large enough. For the dissipative part Z we have, under 
the same conditions, the more precise estimate 

(22) \\dPZ(t)\\ L °o <Cmin{l,t- 1 ~i}E f)+1 , 
for another constant C = C(Ep +a ) as previously. 
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Theorem 2.4. Under the assumptions of Theorem ] 2. SI we have the following decay 
estimates for the time derivatives of Z: 

(23) ||^Z t (t)|| L ~ ^CminU,*- 1 -^}^; 

(24) \\d^Z t (t)\\ L ^ < Cmm{l,t-§-%}E 0+3 ; 

(25) \\d£Ztt(t)\\ L ao <Cmin{l,t-|-t}^ +3 ; 
with C = C(Ep +a ) for a large enough. 

Proof of Theorem (|2.4p . We have only to prove inequality (|23|) , since the other ones 
are proved in First we have that 

(26) d tt u = Al^u + A n A 22 d xx Z - A 12 d tx Z, 

and so, using inequalities ([2"2"]l . (|23|) and (J2U), we have (|2"5"|) for the conservative 
part. Now, using the second part of system dU), we have 

(27) d u Z = -A 21 d tx u - A 22 d tx Z + Q'(u)d t u - d t Z, 

which yields the proof, using that the term Q(u) is at least quadratic in u. □ 

Remark 2.5. Consider the special case of system (0| when An = 0. We can 
improve the decay of the time derivative of the unknowns. Actually 

= ||Ai 2 0xZ|U« < Cminil^-i}^- 

Then, we have 



(28) ||$tu(t)||ioo = \\A 12 d tx Z\\ L - < Cmmil^-^Ei, 

(29) \\d tx u(t)h°° = \\A 12 d xx Z\\ L cc < Cmmi^tr^Ei. 
Finally, for the second derivative of the dissipative part, we have 

(30) d tt Z = -A 21 d tx u-A 22 d tx Z + Q'(u)d t u-d t Z. 



Therefore, we can apply the Duhamel's formula to obtain 

d t Z = e-'dtZo - J e-(*- s > (A 21 d tx u{s) + A 22 d tx Z{s) + Q'(u)d t u(sfj ds. 

Now, since the function Q is quadratic in u and using the previous estimates, we 
find that 

\\d t Z\\ L co <c( y e.- t + J e- (t - s) min{l,s- 2 }(is^ <Cmin{l,^ 2 }. 

By plugging this last inequality in pop we obtain 

(31) \\d t tZ\\ L °° + \\d t Z\\ L °° < Cmin{l,t~ 2 }. 

3. The Numerical Approximation 

In this section we first introduce general finite difference approximations for 
system (fTJ). Then, we compute the local truncation error of these schemes and we 
discuss its decays properties. The main result is given in Theorem l3.1[ where a class 
of TAHO schemes is fully characterized. First, we approximate the differential part 
following the direction of the characteristic velocities, so we study the methods for 
the system in diagonal form (JTJ) . 

We denote by / = (J 1 ,...,/™) the exact solution. Let Ax the uniform mesh- 
length and Xj = j Ax the spatial grid points for all j E Z. The time levels i„, with 
to = 0, are also spaced uniformly with mesh-length At = t n+ i — t n for n £ N. We 
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denote by p the CFL ratio p = At/ Ax, which is taken constant through all the 
paper. 

We consider the Cauchy problem for system ^ possibly subjected to some 
stability conditions. The initial data f° is supposed to be smooth and approximated 
by its node values. The approximate solution (fj, n7 ---,fj n n ) T , f] n G M. mi , i = 
1, ...,m, for jgZ and n G N, is given by 

fj,n+l ~ fj,n , i ^ i x ^ r2 ri 

At 2 Ax /j '- 1 ' n; 2Ax xJj ' n 

(32) = ]T (Si(«i+i,n)-^ +Il „), 

2 = -l,0,l 

/£o = fe), J'eZ, 

where $xfj,n = (fj+l,n — 2/_,. n + /j+i,„), for all i = 1, m. The artificial diffusion 
terms qi are diagonal matrices in R" ilXmi . The source term approximation is 
defined, for I = — 1 , 0, 1, by the diagonal matrices j3\ G W niXmi and by the vectors 
of functions G M mi . 

We assume the scheme (pJ2j) is consistent with system (fTJ), i.e, for all i = l,...,m 

(33) 

S l _!(u) + S' (u) + %{(u) = Mi(u) + AxC, (it), 

where C % = diag(c\, ...,c l mi ) G R m i xm i and Cj(u) are mi functions to be defined. 

By applying the change of variables the scheme applies to the system in the 
conscrvativc-dissipative form and it reads 



Zj,n+i - Zj.n A . 

At + 2Ai {Zj+1 ' n ~ Zj -^ n) ~ 2Ax KZ] n 
(34) 

= ^ ( ^l(uj+l : n) - hZj+l,n) j 
Z= — 1,0,1 

where Q = diag(qi, q m ), Q = DQD^ 1 and for I = — 1, 0, 1, 
k = diag(fil' 1 ,...,0l' k ,...,0P' 1 ,...,0^' k ) G M*™**™, bi = DbiD^ 1 G R km * km , 

SK^) = (^ 1 (n),..,^ fc (n),..,Sr' 1 W ! ...,3 m . fe (n)) T GK fcro ! 
23; = G K fcm . 

By separating the system with respect to the two variables u and Z ', we get 
(35) 

Uj>+1 - Uj, n . Ml , s . M 2 ( ~ ? \ 

2Ax ( - Uj+1 ' n ~ u ^~ 1 ' n ' + 2Ax. \ j+1 ' n ~ ^- 1 ' n ) 



At 

Qll c2 



z=-i.o,i 

Zj,n+l - Zj,n Ml , . ^22 £ 

At + 2A^ tUj ' +1 ' n " Uj -^ n > + 2A^ r i+1 ' n ~ Zj - X ' 

Q21 r2 - 



2Ax 
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where, for I = —1,0,1, vector T$i{u) = (3j(u), 23;(it)) T , follows the notation given 
in an d the block decomposition of matrices bi and Q follows the definition 
given in (fl2l) . 

3.1. Decay properties of the local truncation error. In this section we focus 
on the local truncation error for the general scheme By applying the time 

decay properties given in Section l2~2l we will show how it is possible to build up 
numerical schemes which are more accurate for large times. 
Set, for i = 1, m, 

C = diag{C l ), C = DCD-\ Q{u) = {Q t {u)) T , 7* = (# - 

(36) 

<3 = (51-6-0, r = (S!(u)-s_i(u)), f = (Si(«)-s_i(«)). 

For i = 1, .., m, we have 

%{u j+Un ) = Sj(uj,„) + I Axfyiuj^dvuix^tn) + 0{Ax 2 ), 

where 23 J £ B™" is the Jacobian matrix of 23;. Using the Taylor expansion 
and the consistency property (|3"3"|) . the local truncation error for the scheme (|3"5|) 
becomes 



(37) 

and 

(38) 

where 
(39) 

TO 



T u = -^-C«u - Ax-— d^u - Ax—^-d xx Z 



-Aa 



T u (u) + Tfdzu + S#Z + Sfd x Z + 0{At 2 + Ax 2 ), 



7, 



v a Q21 a Q22 Q ^ 

— d tt z - Ai— d xl u - Ax—^-d xx Z 



+Ax 



T§ (u) + Tf 9 x u + SgZ + 5f a r Z + 0(At 2 + Ax 2 ) 



iJiie^uJ+DiaCCuJ-Cnu) 6 



= -f i + G u e 



t)hxh 



T *(u) = -(DaiC^u) + Z? 22 C(u) - C 2 iu) £ 
SA l = Cia G R fcxm2 , 



C22 € 



Tf = (-r +G2O g R m2Xfc , 

ST = G12 G R fcx ™ 2 , 
= G22 G 



xm 2 



where f is the TO2 X k jacobian matrix of f . Clearly, the scheme is at least 
consistent, which means it is formally of order 0(Ax + At). Taking Ti = 0, 7* = 0, 
C l = and C,; = 0, for i = 1 =,..., m, we get the standard upwind scheme with the 
pointwise approximation for the source term and the local truncation error is just 
given by 

7 u (x, t) = ^d tt u - AxQ±d xx u - Ax^d xx Z + 0(Ax 2 + At 2 ), 



2 

At 



2 
Q21 



Q 



2 9 tt Z - Ax^d xx u - Ax^d xx Z + 0{Ax z + At'). 



7 z (x,t) 

Using the time decay estimates in Theorems 12.31 and 12.41 we obtain for a general 
approximation the following estimates for the local truncation error, as t — > +00, 

7 u {x,t) = 0(Ax t~ 3 / 2 ) + 0{At t~ 3 ' 2 ), 7 z (x,t) = 0(Ax t~ 3 / 2 ) + 0{At t~ 3 ' 2 ). 
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Starting from the general scheme we would like to improve the decay 

property of this local truncation error to build up more accurate numerical schemes. 
The main idea is to chose the free parameters of the scheme to delete the terms 
that decay more slowly in ((37|) -(f38 |) . i.e. the terms which decays as i -3 / 2 . 

Let gi = diag(7 (i _ 1 ) mi+1 , ...,j imi ) for i = 1, ...,m and G = diag(gi, ...,g m ). 

Theorem 3.1 (Local Truncation Error). Let At/ Ax = p be fixed and let H = 
diag{h\, . . . ,h m ) be the block diagonal matrix given in ((5J). Recall that, by 
An = ^ZiLi 7 an d se t P = Sl=i ^i^h ■ Assume An ^ and that the 

following condition holds: 

(40) the matrix {Xilk — An) is invertible for i = 1, m. 

// we make the following choice for the coefficients of the scheme (132[) , 

(41) C = -|/ fcro , S = CM{u) = -|m(u), 



(42) 9i = - QftV 1 - f K 1 (P - (XJ k - An) 2 )j {XA - A n )- X hi 



and 



(43) I» = 9i Ml(u) + P - (V 1 - Ml(u))A n + A,M?(u) - hT l £ ^'(u) , 



both for i = 1, ...,m, then the local truncation error of the scheme (|32p decays as 
(44) 

T„(x, t) = O (Ax t- 2 ) + O (Ax 2 £- 3 / 2 ) , 7 z (x, t) = O (Ax t~ 2 ) + O (Ax 2 i~ 3 / 2 ) . 
Proof. Set 

To = -Q(u), 

(45) fi = -(A 22 (Q'(u)) + (Q'(u))A n - A 21 ), f 2 = A 21 A U + A 22 A 21 , 
5i = 2A 22 - (Q'(u))A 12 , S 2 = A 21 A 12 + A 2 2 . 



Replacing these expressions in ([371) - (l38l) . and using the structure of the system, 

yields 

(46) 

T u = Ax 



T, 



y(u) + (l? + ^(Q'(u))) d x u 



(|a X2 + S?) d tx Z + f -% + jAn^aa - S\A 22 



1, = Ax 



|f (u) + T ») + (Tf + p -fi + (SI + |5x)(Q'(u))) 3 x u 



f T 2 - - (5{ + |Si)A 21 ) c^u + (I J ma + Z - (sf + ^Si ) 



P_S 2 -^f-(St + ^)A2i 



Q(Ax 2 t~ 3 / 2 ), 
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Our choice of the coefficients H\ (u), C l and Cj (it) is made just to cancel the 
coefficients of the slowly terms, i.e. 



(47) 
(48) 

(49) 



To = 0, 
|f o (u)+T o » = 0, 

T? + S?Q'(u) = 0, 



P -A 2 



li 



S z = 0, 



S?A 21 = 0, 



(50) Ti + ft t + (SI + ft)Q'(u) = 0, fti-^f- (SI + ft)A 21 = 0. 
Therefore, the local truncation error reduces to 



7„ 



Ax 



|A 12 + Sr ) d tx Z 



(51) 



Ql2 . P 



A 11 A 22 -S^A 22 ) d xx Z + 0(Ax 2 t~ 3 / 2 ), 



7 Z 



Sf + ft) 



^S 2 - % - (Si + ^S!)A 22 ) d xx Z + 0(Ax 2 t^' 2 ) 



and by the estimates given in (|2"T]) - ([2"4"1) the thesis is achieved. 

We need now to show that system (|4"7l) - ([50|) has always a solution given by the 
relations (|4"Tj) - (|4"3")) . The terms given in (|4"Tj) for the C and C coefficients are obtained 
from 07} an d (03). Indeed, we get 



-DiiG^u) - £>i 2 e(u) + Cnu = 



C 



12 



0. 



-Q(u) + (-^iC 1 (u) - D 22 G(u) + Cam) = 



c 2 



or in a more compact form, 



D n D 12 
D 2 i D 22 



e l (u) 
e(u) 



C21 C22 



which can be rewritten as 

-DG + DCD^Z 









Q(u) - I m2 Z 



Q(u) - Z 

Now, we multiply on the left by the matrix D~ l to obtain 



(52) 







Q(u)-Z J 2 



£(-/ + M(t0), 



that gives (|4T1) . 

We shall now focus on the derivation of the relations (|4"2"]l and (|4"3")l . From (|49l) - 
(f50l) , for the matrix of free coefficients G we have to impose 



(53) G 12 A 21 = -^f + ft 2 llt G 22 A 21 = -^f + P -[ 
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Now, using the notation (fT3f for the two diagonal matrices G and Q and by 
applying relations © and ©, we have that for G = DGD^ 1 and Q = DQD^ 1 

G\2 = G\H 1 D 21 + Dr 2 G2H 2 X D 22 , 
G22 = D2\G\H 1 1 D 21 + D 2 2G2H 2 D 22 , 
Q11 = Q\H 1 1 + Di 2 Q2H 2 1 Dj 2 , 
Q21 = D 21 Q 1 H^ 1 + D2 2 Q2H 2 1 Dl 2 . 



Since Dn = Ik, we have that (|53|) becomes 
(54) 

G 2 H 2 - 1 D^ 2 A2 1 J 2 V Q2H 2 1 D{ 2 J + 2 I A 21 A n - A 22 A 21 



where we used relations (|45j) to sort out the term T2 — S1A21 = A21A11 — A 22 A2i. 
Here we used the fact that, neglecting the terms with a faster decay, we can replace 
S\ with 2^22- Actually, since the term S\ in (|46[) multiplies only terms that decay 
faster than i -3 / 2 , it is possible to write 

Si = 2A 22 - Q'{u)A 12 = 2A 22 + 0(Axtr 2 ). 

From ([5"4")) we get 

(55) 

GiH^D&An = -\QiH^ + P -H^ x (A 2 n + D T 21 (A 21 A n - A 22 ^ 2 i)) 

G 2 H 2 - 1 D 2 r 2 A 2 i = -\Q 2 H 2 x Dl 2 + P -H 2 X {D T X2 A\ X + D T 22 (A 21 A n - A 22 A 21 )) , 

Using the specific form of D and relations (fTDj) . by algebraic considerations we 
get 

A\ x + D T 21 {A 21 A U - ^22^21) = P - (Ai/fc - A n ) 2 , 

and 

Df 2 A 2 n + D 22 {A 21 A n - A 22 ^2i) = Df 2 P - Dl 2 A 2 n - K 2 Dj 2 + 2A 2 Df 2 An. 
Then, we get for i = 1, .., m 

(56) gih^iXih - A n ) = -\qJ\ x + ^hj 1 (P - (\J k - An) 2 ) . 

Assuming for i = 1, m, (Ai/fc — An) to be invertible, we obtain relations (|4"2"1) . 

Finally, we need to compute the vector function T(u). From (I50|) we obtain the 
two following relations 

(57) - fi + Gn + G 12 Q'(u) = 0, - f ' + G21 + G 22 Q'(u) = - P - (Tj + SiQ'(u) 



-1 I u \ a r'l„,\ _ n-l I Ik 



M W = D ' Q(u) ) = D ~ 1 Q'(«) 



Since 



equations ([57)) reduce to 

(58) -r'w + Gw'f^-y 1 ' 



2 V T i + W(u) 

where f x + SiQ'(u) = A 22 Q\u) - Q'(u)i4n + A 2X . Since Q'(u) = D 2 iM((u) + 
D22M'(u) and M{(u) + Di2M'(u) = Ik, by using relations (fTU)) , we obtain that the 



12 



D. AREGBA-DRIOLLET, M. BRIANL AND R. NATALINI 



right side of (f5T|) is equal to 
/ 



Hr 1 y(h - H 1 M[{u))A u + A 1 J? 1 M((u) - J2 ^M-(u)J 

/ m 

H 2 X Df 2 A n - H 2 M'(u)A n + A 2 H 2 M'(u) - D T 12 ^ A 4 M/( 



\ \ 1=1 / / 

Therefore, we get for every vector Ti(u) £ Mr, i = 1, m, that 
(59) 

/ m > 

- rJCu) + 9i Ml(u) = - P - (hr 1 - M[{u))Avi + \M[{u) - hr 1 £ X.M^u) 



□ 



Remark 3.2. As previously observed in Remark 12.51 when A\\ = 0, we have, for 
the second-order time derivatives, 

d tt u~t- 2 , duZ^r 2 . 

Therefore, in dealing with (|37p - (|3"5)l , we do not need to delete the second-order time 
derivatives Utt and Z t t, and relations (|4"Tj) - (|4"3"]) reduces to 

(60) C = 0, 6 = 0. 

Besides, for each i = 1, m such that A.; 7^ 0, we can choose 

(61) 9i =fi- = -^qu Ti(u) = -B\ - Sii = 9iMi(u). 
Then, we can select qi = |Aj|,for i = 1, ...,m, and we are free to choose 

(62) /3 = -, T> (u) = 2 ■ 

Therefore, we obtain an upwind scheme for system (JTJ), with the classical Roe 

upwinding approximation for the source term [21) . namely 

(63) 



fj,n+l fj,n , ( fi r i \ I'M , 

i 



A< 2Aa- ^"W 2Ax ^ 

_ Mi^J + 2Afi(u?) + Afi(«J +1 ) S5 n(A,) , 



+ 2 /j,» + sg»(A a ) „ „ 

4 4 Uj— l,n Jj + l,n)i 

From now on, we shall refer to scheme as the ROE scheme. Then, we have 
just proved the following result. 

Proposition 3.3. For An = 0, the local truncation error of the ROE scheme (|63[) 
verifies the time asymptotic estimate f|44[) . 



4. A MONOTONE TlME-AHO SCHEME FOR THE 2x2 CASE 

In this section we apply the general result stated in Theorcm l3.1l for the 2x2 case 
described in the Example (|2.ip . For the artificial viscosity, we assume q\= q 2 = q, 
q e M+. 

Hence, by applying the results given in Theorem 13.11 we get for the scheme 
parameters the following expressions: 

(64) Ci=C 2 = -|, = -|Mi(«), e a (u) = -|M 2 (u) 
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(66) T 1 (u)= 2a M r {u)-—pu, T 2 (u) = — M 2 (u) + -^pu. 

To complete the definition of scheme ([34]) it is still necessary to choose four more 
free parameters, such as 23 ' (■) and j3 ' . For the 2x2 case such parameters can 
be defined by applying the monotonicity conditions. 

Starting from the scheme written in its diagonal form (|32p . the monotonicity 
conditions are given by the following relations, see pQ: 



!B' 1;2i (-)>0 VZ = -1,0,1, 
1-pq + At (V 1;2 , ( U )-/3 1;2 ) >0, 



^ + f + At (81^(1.) -^?)>0, 

"^ + f + At(s' 1;2 »-^ 2 )>0. 

Proposition 4.1 (Monotonicity). Assume a > and A > H-F^u)!^. Under 
the assumptions of Theorem \3.1[ the scheme (|34[) /or i/ie 2x2 case verifies the 
monotonicity conditions (|67[) /or i/ie choices: 



Sj(u) = Afi(u) - |r x (u)| + Ax C 1 ^), S§(u) = M 2 («) - |r a (u)| + Ax C 2 (w). 

(69) $ = 1-71+ Ax ci, $ = I + 72 + Ax c 2 , 

under the CFL conditions 

1 - Ap 1 — Ap N 



(70) At < 

and, for X > 2a 



1 + 7i ' 1-72 



a 2 fl 2A 2 m_ 2A 2 m+ 

(71) Ax < , p < mm — , — — 5-, — — 5 

v ' ~ A + a V A 2a 2 Am_ + a + a 2 _ 2a 2 Xm + + a\a- 

where mi ;2 = min u (M{. 2 (u)) > 0. Otherwise, if a < A < 2a, we get the 
supplementary requirement 

\X-2a\ 

P> 2 A ■ 

a z — Ax a_ 

Proof. From consistency (1331) . we write 

^(u) = \(Mi-, 2 (u) - %^{u) - r 1;2 (u) + Ax Ci ;a («)) , 

®i = J(Mi ;2 («) - ®J ;2 («) + r 1;2 ( u ) + Ax Ci i2 («)) , 

P 1 - 2 ! =l( 1 ~ ti 2 - 7l;2 + Ax c 1;2 ) , ft 1 ^ = i (l - ft 2 + 7l;2 + Ax c 1;2 
The first condition in (|67[) is equivalent to 

(72) M(. 2 - |r' 1;2 1 + Ax C 1;2 > SU 2 > 0, 
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and the third condition in (|67p is equivalent to 

^(±\ + q)-^-(l-plTl2 + &X c+) >0, 



(73) 



£(TA + g)-^(l-$=F7i+Aa: c_) > 0. 



It is natural to assume that the CFL ratio p verifies the standard hyperbolic 
condition 

(74) p<\. 

Then, from relations (|65[) . we have 7 + < and 7~ > and for q > A, we get 
monotonicity by choosing in ([73")) . 

(75) $ = 1 - 7i + Aa; c_, /?§ = 1 + 72 + Ax c+, 

under the limitation required in the second condition in (|67p. 

1 - A/? 1 - Xp" 
1 + 71 ' 1 - 72 



(76) At < n 

To verify the request l|72p. we set first 



(77) S' 1)0 («) = M[ - \T[\ + Ax C[, = M 2 - \V 2 \ + Ax C 2 . 
For 

/ 2A 2 ra_ 2A 2 m + 
p < mm — — T , — — = 

\2a z Am_ + a+al 2a^Am + + a^_a_ 

we get 

r 2 <o, r;>o, 

then conditions (|72|) become 

'-sr<£-*-)Kih* 

that are verified under the following limitations, 

(78) A >max(||F / («)|| oo ,2F / (0)), Ax < 



A + a' 

If we choose 

a < A < 2a, 

we get a limitation from the bottom for p, 

\X-2a\ 
P> 2 ' ' ■ 



□ 
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The monotone Timc-AHO scheme has then the following expression, 
(79) 



</ .7,71+1 J j, 7 



— if 1 - f 1 



(l - ^) ^Kn) + ^M 1 {u i+hn ) ((1 - 7 -)/l„ + l-f} + l, n ) 



f 2 - f 2 



At 



= ^M&i-hn) + (l - ^T^) M 2 ( Uj . n ) - (| 7 + |/jU,„ + (1 - l7 + l)/|„) 

+f§( Uj , n - Uj ^ n ) + £§2 (f 2 n - M 2 ( Uj , n )) . 

Notice that, the approximation of the source terms involve only the values of the 
solutions on the " upwinding- nodes" , i.e. (xj,Xj+i) and (xj-\,Xj) respectively for 
the first and the second equations. 

The scheme may be then considered as an extension of the known approximation 
with the upwinding of the source term (|63|) described in Remark 13.21 Let us stress 
on the fact that, when F'(Q) — 0, the local truncation error of the ROE scheme 
((63|) verifies the decay properties obtained by the Time-AHO schemes in (l44l) . 

5. A monotone Time-AHO scheme for the 3x3 case 

In this section we compute a TAHO approximation for the 3x3 case of example 
(|2.2[) . In particular, we study the case where F(u) = a(u — u 2 ), with a > 0. Then 
we have F'(0) = a > 0. 

We recall that a = \a\ + j3X e]a, A[ and 

(80) X-a+\a\>\F'(u)\. 



According to Proposition 13. 11 we have An = a, and P = Xa, and the coefficients 
hi are 

2A A , 2A 
h x = , h- 2 = , h 3 = — ■ — . 

a — a A — a a + a 

We choose to take q\ = q% = A. q 2 = 0. Consequently we find 
(81) 

A(l — pa) p(X + a) p(Xa — a 2 ) A(l — pa) p(X — a) 

9i = oM , x H o ' 92 = ^ , 93 = 



2(A + a) : 
We have 

r x (u) 

(82) \ r 2 (u) 

r 3 («) 

For 1 < i < 3 we have 
1 



2a 



2(X-a) 



A(l — pa) , r . . p(a~a), 

; -Mi(u) + ^ - ' {au - F(u)), 



2(X + a) 
pXa 



2a 



M 2 {u) + 



4A 

p(X - a) 



2X 



(au — F(u)), 



A(l — pa) , . , . p(a + a) . . . . 

v P ; -M 3 (m)+ PV - ' (au~F(u)). 



2(A - a) 



J ±i 



o > 



2 



±i 



4A 



A/ 



)±r 4 («)-!B*(«) 
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so it remains to determine the /3q and H^u). As for the 2x2 case we use 
monotonicity criteria for those choices. From now on, we use the fact that a > 
and we consider only u € [0, 1], as this will be satisfied in our numerical test. We 
have then 

(83) -a<F'(u)<a, 

a , , . a — a 

(84) < 1 - - < M' 2 (u) < 1 - —j- < 1, 

(85) <^<M/K><^<1, i = l,i = 3. 

Hereafter we study the monotonicity conditions for each of the three equations 
obtained by applying the general numerical scheme (|32[) to this particular 3x3 
case. We denote 

fi = (Ax + 2a) (X + a) - Xa, 

and 

A / a — a\ Xa ( a — a\ a(a + a) 

" 1 = 2(X—a){ 1 -—)> " 2 = -—{ 1 -—) + ^^ + X - a - 

Proposition 5.1. (Monotonicity). Assume a > 0, A > 2a. We suppose that the 
following conditions are satisfied: 

(86) if v 2 > then Ax < min 



(87) if v 2 < then Ax < min 

(88) p < min 



A + a> + ^;' 
Aa A 



A + a ' v\ ) ' 

11 2a 2a(X-a) 



a ' A + Ax ' aAx + Xa ' a(X - a) Ax + a(X 2 + A - 2a 2 ) 

X + 2a 



(89) if p > then p < 

(90) At < 



1 



1 + max \T' 2 (u) — g 2 \ 
ue[o,i] 



Under the assumptions of Theorem \3.1[ scheme (|34[) for the considered 3x3 case 
is monotone for the choices: 



At „ .At „, At 

a ±7 r- u ~ 2 ' Ji! / " u ~ 2^ ^ ' 



(91) ^ = 1- — -ffi, /3 2 -l- — +.g 2 , /3 3 = 1 



(92) S5(«) = Mx(«)(l - ^) - ri(u), S§(«) = M 3 (u)(l - ^) + r 3 («), 
and 23q is a continuous function such that 

(93) - (1 - M£(u))(l - ^) - \T' 2 (u) g 2 \ = (Sg)'(u) - ft. 
Such a function exists. 

Proof. We do not detail the proof for equations 1 and 3, as it follows the lines of 
the 2x2 case. Actually, under the above assumptions we have 

gi >o, ri(u)>o, 33 <o, r' 3 ( u )<o. 

Then we obtain a monotone TAHO scheme on those equations. 
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The treatment of the second equation is quite different. We first note that g 2 < 
but the sign of T' 2 (u) does not depend on the parameters of the discretization. The 
monotonicity conditions are: 

(94) CB?)'(u)>0, 1 = -1,0,1, 

(95) -0 2 ±1 + (03 2 ±1 )» > 0, 

(96) 1 - At(f3 2 - (Sg)'(«)) > 0. 
The inequality (1951) can be written as 

(97) - (1 - M 2 ( U ))(1 - ^) - |I" a («) - ^1 > CBg)'(«) - /3 2 , 

which is implied by equality (|93|l . In the case under consideration, it is 
straightforward to determine the sign of T' 2 (u) — g 2 with respect to u. In our 
numerical experiment for example, we have uq €]0,0.5[ and u\ s]0.5, 1[ such that 

T' 2 (u)-g 2 >0 in [0,u [U]ui,l], 

r 2( w ) - 92 < in ]«o,«i[- 

We can therefore construct a continuous function 25q such that 

CB 2 y(u)^M 2 (u)(l-^)-T' 2 (u) + 2g 2 when r' 2 ( U ) - ff2 > 0, 
(Sg)'(u)=M^(«)(l-^)+r' a (u) when r 2 (u) - ff2 < 0. 



We set 
Therefore 



P-l = -32, 0Q=l-Y+to> /? 1 =0 - 



-(1 - M£(u))(l - ^) - |I» -g 2 \ = (Sg)'(u) - /3 2 < 0. 



Consequently we have (|95p. (f9"6")l by (|90p. and (|94p for I = ±1. It remains to satisfy 

(98) (S 2 )» > 0. 

First case: r 2 (u) — g 2 > 0. 
By (P |) -([54 ll we have 

A — a , . 

> -jT-(l-<7p). 



One can prove that <r > 0. We obtain ([98]) by 
Second case: T 2 (u) — g 2 < 0. 

(BS)'(u) > M ^)(l"^-^)>0 
by ((88|). □ 

6. The linear case 

6.1. Numerical schemes. Here we want to apply the argumentation of the above 
sections to the linear case, first considered in work [T|. We shall focus on problem 
(|T7|) for F(u) = au. We set a = jia and (3 = \x and we study the following problem 



(99) 



u t + au x + z x = 0, 
zt + u x - az x = -j3z. 
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The numerical approximation given in (|34[) , here becomes 

T n+l _ 

U -J- 4- 

2Ax 



( 100 ) At 2Ax v J+1 2A.t 

= B.xC/J 1 .! + B lH l + Sit^ +1 , 

where U = (u, z) T , for ±£ = ±yl + a 2 be the eigenvalues of matrix A, Q = 
diag({;,£) and ®_i,o,i = (/^ ' ' )j,j=i,2 are the matrix of constant coefficients for 
the source term approximation. 

First of all we shall analyze the decay properties of the local truncation error for 
the numerical approximations described in pQ. 

By Taylor expansion, the local truncation error for the numerical approximation 
(|100P is given by 
(101) 



T 



7, 



¥u tt - Ax 



U X x + Wli ~ Pn)Ux + (012 ~ Pl2~) Z x + C 11 U + C 1? Z 



+0(Ax 2 + At 2 ), 



^rzu - Ax 



Zxx + (021 - 02l) U x + (022 ~ 022 1 ) Z x + ^2\U + C 22 Z 



+0(Ax 2 + At 2 ), 

where the (cy)i,j=i,2 constants were defined in [1]. Let At/ Ax = p be fixed. By 
relations 

u tt = a 2 u xx - (z t - az x ) x , z tt = £?z xx - 0(z t + az x ), 
we get, for some popular schemes, the following expansions. 



(UP) for the source term point wise approximation, we have 
(102) Si-S_ x = 0, C = 0, 



(103) 



T 



7, 



^[(pa 2 -t)u xx - p(z t -az x )x\+0(Ax 2 ) ~ Axt~ 3 / 2 , 



^ - l)z xx - P0(zt+az x )} + 0(Ax 2 ) ~ Axt- 3 / 2 . 



(ROE) for the upwinding of the source term, we have 

( 1 



(104) 



C = 0, 



(105) 



2£ \ -a /' 

( pa 2_ i^2l )Uxx + a_ p )( Zt -az x ) J +0(Ax 2 ) ~ Axt- 3 / 2 , 



ry . Ax 

z 2 



Z(pt - l)zxx - 0(pz t + a(p - ±)zx)\ + 0(Ax 2 ) ~ Ax i" 3 / 2 . 
(AH02p) for the Asymptotic High Order scheme given in [T], we have 

2 t ( 



(106) 
(107) 



m m ft ( 1 



1 



7 U = [pa 2 u xx + (€ - p) (2* - a^)J + 0(Ax 2 ) - Ax i" 3 / 2 , 



A.r 

2 



Upti-l)zxx- 0(p + O(zt+az x )] + O(Ax 2 ) ~ Axt" 3 / 2 . 
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Let us now go back to the Time-AHO schemes. According to the discussion of 
section 13.11 conditions given in propositions I3.1M.11 here give rise to the following 
choice, 



For the local truncation error it then holds, 
(109) 



7 U = K(l - P 0{-Zxt + az xx )] + 0{Ax 2 ) ~ Axt~ 



7 Z = -^[Z(l-p0z xx } + 0(Ax 2 ) ~ Air 2 . 

Remark 6.1. As for the non-linear case, for a = 0, the Time-AHO scheme reduces 
to the Roe approximation, which in that case decays like in (|109p . 

6.2. Modified equation. Here we wish to better understand the qualitative 
behavior of the numerical methods described in the previous section by applying 
the modified equation method (see for instance |18j). From Taylor expansion stated 
in (|101j) . for all numerical schemes described in Section H3 we get the following 
modified system 
(110) 



u t + au x + z x = - ^rUtt + Ax 



z t +u x -az x +f3z = — -^z tt + Ax 



^u xx + juUx + l\iz x + c n u + c 12 z 

\z xx + 72lW x + J22*x + C 2 lU + C22Z 



Since we are interested in long time simulations, we apply the decay rates results 
given in the analytical Section |2~21 to the derivatives terms of problem (|110[) . As 
time increases, we can then take into account as modified equation the following 
asymptotic modified system, 
(111) 



u t + au x + z x = -^u tt + A 



§U XX + 7111*3 + 7122a + CllW + C i2 Z 



2 

[ U x + j3z = Ax[j2lU x + C2\U + C22Z] . 

For all schemes described in Section [5J we get in the second equation 
[721 + C21U + C22z\ = and then, for all of them the asymptotic modified problem 
becomes 

(112) 



u t + au x = U + AxD)u xx + 0(Ax 2 ), 
z= -±u x + 0(Ax 2 ), 



where the constant D depends on the selected scheme. 

Specifically, for UP, ROE and AH02p we obtain a perturbation of order O(Ax) on 
the diffusion term with 

m 1 t 2 t \ m 1 ( 2 C 2 - A ^ pa 2 
^>UP^--^{pa -£), T) RO B = 1 — )' Ti aho2 P = - — ; 

while, for TAHO we get 

'Dtaho = 0. 

The TAHO approximation is then second order accurate on the asymptotic 
diffusion (Chapman-Enskog) limit 
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7. Numerical tests 

In this section we show how, for large time simulations, Time-AHO schemes give 
better numerical results than standard approximations both for linear and non- 
linear cases. 

From now on we shall call STD the following standard pointwise upwind 
approximation, for i = l,...,m 

(114) [ln±l - [in , A ' / fi _ fi \ _ x2 p _ M J . \ _ fi 

For all tests, we focus our attention on the numerical error as a function of time: 
denoting (u,Z) the conservative-dissipative variables, we plot the errors e u (t) = 
\\(u H -U h )(t)\\ L ~, e z (t) = \\(Z H -Z h )(t)\\ L °° as the time t = nAt increases, where 
(u H , Z H ) is the reference solution obtained by the ROE scheme with Ax = O(10 -4 ). 

For all schemes, we fix the steps ratio p to verify all the CFL conditions; Since 
all schemes are of first order approximation, to emphasize the good behaviour of 
TAHO compared to the others schemes, we compute the numerical solutions U h by 
using a quite big grid step Ax = O(10~ 1 ). 

We then plot the different approximations of functions u and Z at final time 
T = 450, focusing on the point of maximum value of the solution to highlight the 
differences of the approximations. Near it, we show the most interesting plot of the 
l°° errors as a function of time. 

Then, given different numerical approximations U , we look for constant C u , j u , 
C z , 7 Z which best fit the equality 
(115) 

e u (t) = \\(u H ~ U h )(t)\\L- - a*" 7 ", e z (t) = \\(Z H - Z h ){t)\\ L ~ = G z t~^. 

Given N data points (ij, e(ij))j = i i jv ) we shall define 7 and C as the solution of the 
following least squares problem, 

JV 

min^|ln( e (t,))-ln(C^)| 2 . 

' 7 i=l 

All numerical results we present show that for standard approximations, such as 
upwind (|114|) and ROE |63l) . the absolute error e(t), for a fixed space step, decays 
as 

e u (t) = Oir 1 / 2 ), e z (t) = Oit- 1 ), 
while for the TAHO scheme, it improves of t^ 1 / 2 on the previous schemes. 

7.1. Results for the linear case. Let us consider for system ([M]) . the constant 
equilibrium state u = 1 and z = 0. As in our previous work pQ, we fix fi = 5 and 
we consider a small compactly supported perturbation of this constant solution as 
initial data. 

(116) m = X[-i,i] (-x 2 + 2) , z = Xl-1,1] {-x 2 + l) • 

We then compare the Time-AHO scheme (|108|) with the following schemes: the 
AH02p scheme p06[) , the standard first-order point wise upwind scheme (|102p and 
the ROE scheme (fTMj) . 

7.1.1. Test case with a = 1. As expected by our asymptotic analysis, the numerical 
results show a better performance of the TAHO scheme. In Figure[T]-(a)-(b) we plot 
a zoom on the solutions u and z respectively, obtained by applying the different 
numerical schemes at final time T. The solution given by TAHO follows much 
better than the other the benchmark curve. 

Moreover, always in Figure [TJ the two plots (c) and (d) show the time evolution 
of the l°° errors e u (t) and e z (t) defined in (| 1 1 5[) for all schemes considered; They 
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show how for the TAHO scheme, as time increases, both errors decay more quickly 
than the other. This result is also confirmed by Table [TJ where the values of 7 
and C are computed. Looking at the different values of 7, it is clear that for the 
TAHO approximation the decay velocity of the absolute error improves of i -1 / 2 on 
the previous schemes. 

To stress on the good behavior of the TAHO scheme, in Figure [21 we plot the 
solution u obtained by different numerical approximations with decreasing space 
step Ax. All schemes considered are of first order approximation, but looking at 
the numerical curves, it is clear how for large step the TAHO solution follows better 
the benchmark line than the other. 
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Figure 1. Linear Case Test, section 17.1.11 (a)-(b) zoom on the 
solutions u and z respectively obtained by the different schemes at 
final time T. The plot show a better performance of TAHO scheme. 
(c)-(d) time evolution of the l°° errors e u (t) and e z (t) defined in 
(|115[) for the different schemes. As expected by our asymptotic 
analysis, for the TAHO scheme the absolute errors e UyZ {t) decay 
faster as the time increases. 



7.1.2. Test case with a = 0. As previously observed in Remark 16.11 for a = the 
ROE scheme (| 1 04[) corresponds to our TAHO approximation. 

For this particular case we can compare the TAHO/ROE scheme with the well- 
balanced approximation proposed by Gosse and Toscani in jlOj . which however is 
defined only in the case a = 0. From now on we shall refer to this scheme as 
WB-GT. 

Figure [3] shows the performances of TAHO/ROE, STD and WB-GT scheme for 
problem 

u t + z x = 0, 
z t + u x = -£lz, 



(117) 
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scheme 


Cu 


7u 


c z 


7z 


STD 


0.192993 


0.510798 


0.100412 


1.058824 


ROE 


0.143847 


0.573922 


0.073287 


1.112843 


AH02p 


0.104897 


0.547166 


0.072874 


1.143070 


TAHO 


0.289768 


0.961310 


0.141281 


1.432194 



7 and C for e u (t) ~ C u t" lu and e z (t) = C z t~ lz defined in 
(|115[) . For standard approximations, such as STD and ROE 
, the numerical results show that the absolute error decays as 
e„(i) = 0(^ 1/2 ) and e z {t) = O^" 1 ); while, for the TAHO scheme, 
it improves of t^ 1 / 2 on the previous schemes. 
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Conservative solution for Dx = 0.0125 
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Figure 2. Linear Case Test, section [7.1. II Zoom on the solutions 
U obtained by the different schemes at final time T by applying 
decreasing values of Ax: (a) Ax = 0.1, (b) Ax = 0.05, (c) Ax = 
0.025 and (d) Ax = 0.0125. The numerical solution obtained by 
TAHO scheme follows better than the others the benchmark curve 
already with a quite big space step. 



with P = 5, at final time T = 450. 

We observe that both TAHO /ROE and WB-GT give better performances than 
STD. The two numerical solutions obtained by TAHO/ROE and WB-GT showed 
respectively in Figure 0-(a)-(b) are overlapping and the numerical errors in Figure 
|3]-(c)-(d) have a similar trend. 

Indeed, according to the discussion of section I3~T1 for the local truncation error of 
the WB-GT approximation it can be shown that 

T„ - Ax t~ 2 , 7 Z - Ax t- 3/2 . 
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Then, we observe that the WB-GT scheme is Time-AHO with respect to the 
conservative variable u, while for the dissipative one it is not. The good behavior 
of WB-GT scheme may be confirmed by analyzing its order of convergence with 
respect the asymptotic modified problem (|112[) . As done in section HOI when t goes 
to +00, it is possible to show that the WB-GT scheme is second order accurate 
with respect the asymptotic modified problem 

L j3Ax 2 
at ~ ~^ xx + 4(1 + /3 Ax/2) V 

z 



namely it is of second order with respect the Chapman-Enskog limit. Notice that, 
for a 7^ 0, it is possible to consider other Well Balanced scheme as in [JJ, but they 
are not THAO and their asymptotic performances are not of the same order (not 
shown) . 
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Figure 3. Linear Case Test for a = 0, section [7X^1 (a)-(b) zoom 
on the solutions u and z respectively obtained by the different 
schemes at final time T with Ax = 0.08. The plot show a better 
performance of TAHO/ROE and WB-GT scheme, with respect 
the STD approximation. As described in section 17.1.21 WB-GT 
scheme is Time-AHO with respect the conservative variables and, 
as TAHO/ROE, it is of order 0(Ax 2 ) with respect the parabolic 
asymptotic state, (c)-(d) time evolution of the l°° errors e u (t) and 
e z (t) defined in (|115JI for the different schemes. As expected by our 
asymptotic analysis, for the TAHO/ROE and WB-GT schemes, 
the absolute errors e u<z {t) decay faster as the time increases. 



7.2. Results for the non-linear 2x2 test case. We fix q = A and we compare 
for the 2 x 2 case the scheme TAHO (7J|) with ROE flM} and STD ([Hi]) . 
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We shall consider the two different cases F'(0) 7^ and F'(Q) = and we select as 
initial datum the function 

(119) u = X[ -i,i] (-x 2 + 1) , Z = ~F(u (x)). 

7.2.1. The case F'(0) ^ 0. Here we fix 

F(u) =a(u- u 2 ) 

and we compare our TAHO approximation (|T9"| with STD and ROE schemes, 
defined in (|114[) and (j6"3")l respectively. 

The numerical results show a better performance of the TAHO scheme. In Figure 
H]-(a)-(b) we plot a zoom on the solutions u and Z respectively, obtained by the 
different schemes at final time T. The solution given by applying the TAHO scheme 
follows much better than the other the benchmark curve. We stress on that the 
numerical solutions are computed with quite big step Ax = 0.1. 
Moreover, always in Figure |H the two plots (c) and (d) show the time evolution 
of the l°° errors e u (t) and e z (t) defined in (|115|) for all schemes considered; They 
show how for the TAHO scheme both errors decay as time increases more quickly 
than others. This result is also confirmed by Table [2J where the values of 7 and C 
are computed. Looking at the different values of 7, it is clear that for the TAHO 
approximation the decay velocity of the absolute error improves of t~ x l 2 on the 
previous schemes. 



scheme 




7« 


c z 


lz 




STD 


0.013797 


0.374708 


0.010744 


0.341554 




ROE 


0.004874 


0.333634 


0.007850 


0.439996 




TAHO 


0.111380 


1.151517 


0.495480 


1.451030 





Evaluation of constants 7 and C for e u (t) = C u t~ lu and e z {t) = 
C z t~ lz defined in (|115|) . For standard approximation STD and 
ROE, the absolute error decays as e u . 2 (i) ~ Oft™ 1 / 2 ); while, for 
the TAHO scheme it improves of tr 1 / 2 . 



7.3. Results for the 3x3 system. As initial data, we take the smooth function 
Mo defined by 

u (x) = xi-1,1] cxp - Y~^) ■ 

Then we set fo(x) = M{uq{x)). We know that in this case one has u(x, t) € [0, 1] 
for all (x,t) e R x [0,+oo[. We choose a = 1, A = 2.1, /3 = (a - a)/X = 0.1. The 

discretization parameters are Ax = 0.1, p = — , which satisfy all the monotonicity 
requirements, see proposition l5.ll 

The numerical results show a better performance of the TAHO scheme. In 
Figures [5j |6]-(a)-(b), 0-(a)-(b), we plot the solutions u, Z\ and Z% respectively, 
obtained by the the STD, ROE and TAHO schemes at final time T, as well as 
the exact (reference) solution. The solution given by applying the TAHO scheme 
follows much better than the other the benchmark curve. We stress on that the 
numerical solutions are computed with quite big step Ax = 0.1. 
Then in Figure [8]- (a)- (b), we plot the time evolution of the l°° errors e u (t) and 
e z (t). They show how for the TAHO scheme both errors decay as time increases 
more quickly than other. This result is also confirmed by Tabled where the values 
of 7 and C are computed. Looking at the different values of 7, it is clear that for 
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Figure 4. Non-Linear Case Test with F'(0) 7^ 0, see section [7XT1 
(a)-(b) Zoom on the solutions u and Z respectively obtained by 
the different schemes at final time T. The plot show that TAHO 
scheme gives better results than others with a quite big step Ax = 
0.1. (c)-(d) Time evolution of the l°° errors e u {t) and e z (t) defined 
in (|115[) for the different schemes. As expected by our asymptotic 
analysis, for the TAHO scheme the absolute errors e u ^ z (t) decay 
faster as the time increases. This result is confirmed in Tabic 
[21 where we compute the decay parameters 7 of absolute errors 
previously plotted. 



the TAHO approximation the decay velocity of the absolute error improves of t 1 / 2 
on the previous schemes. 



scheme 


C u 


In 


c z 




STD 


0.0052 


0.54 


0.0064 


1.1 


ROE 


0.0027 


0.66 


0.0036 


1.2 


TAHO 


0.006 


1 


0.012 


1.62 



of constants 7 and C for e u (t) = C u t~ lu and e z (t) = C z t u 
defined in (| 1 1 5[) . For STD and ROE approximations, the numerical 
results show that the absolute error decays as e u (t) = <3(t -1 / 2 ) and 
e z {t) = Oit- 1 ); while, for TAHO's it improves of i" 1 / 2 . 
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Figure 5. The 3x3 system with F'(0) ^ 0. Left: u component 
at final time T. Right: detail. The reference and the TAHO 
computed solutions cannot be distinguished. 
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FIGURE 7. The 3x3 system with F'(0) ^ 0. Left: Z 2 component 
at final time T. Right: detail. The reference and the TAHO 
computed solutions cannot be distinguished. 
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